Dynamics of Random Graphs with Bounded Degrees 
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We investigate the dynamic formation of regular random graphs. In our model, we pick a pair 
of nodes at random and connect them with a link if both of their degrees are smaller than d. 
Starting with a set of isolated nodes, we repeat this linking step until a regular random graph, 
where all nodes have degree d, forms. We view this process as a multivariate aggregation process, 
and formally solve the evolution equations using the Hamilton-Jacobi formalism. We calculate the 
nontrivial percolation thresholds for the emergence of the giant component when d > 3. Also, 
we estimate the number of steps until the giant component spans the entire system and the total 
number of steps until the regular random graph forms. These quantities are non self-averaging, 
namely, they fluctuate from realization to realization even in the thermodynamic limit. 
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I. INTRODUCTION 

A random graph is a set of nodes that are connected 
by random links fll-[3| • When the number of links exceeds 
a certain threshold, a giant component with a finite frac- 
tion of all nodes emerges [1, Q . Therefore, random graphs 
are equivalent to a mean-field percolation process 0, Wi ■ 
Random graphs underlie many natural phenomena from 
polymerization to the spread of infectious diseases 

12 - 14 1 , and they are also used to model social networks 
15 - 17 1 as well as complexity of algorithms [l^ . 
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The classical random graph has no restrictions on the 
degree of a node. In many situations, the number of con- 
nections is limited, and specifically, the degree of each 
node must be bounded. Examples include the percolat- 
ing network of contacts in a bead pack |19l|. c omputer 
networks [2^ , and communication networks [2ll | . 

We modify the canonical model for an evolving random 
graph 0-11, [23 - [23 by imposing a strict bound on the 



degree of a node 



26l |. In an evolving random graph, 
two nodes are picked at random, and regardless of their 
degrees, the nodes are connected by a link. At the end of 
this process, a complete graph forms. In this study, we 
focus on a constrained process where the two nodes are 
linked only if both degrees are smaller than d. Starting 
with a set of N isolated nodes, the final state is a regular 
random graph [27, 28], where all nodes have degree d 
(Fig.[T]). The evolving graph includes two types of nodes, 
active nodes with degree smaller than the maximum, and 
inactive nodes with the maximal degree. We obtain the 
density of active nodes from the full degree distribution, 
which is shown to be a truncated Poisson distribution. 

In general, the graph is a set of disjoint connected com- 
ponents [29|. In each component, any two nodes are con- 
nected by a path. As in classical random graphs, our 
graphs could be in one of two phases — a non-percolating 
phase where all components are finite, and a percolat- 
ing phase where a single giant component with a macro- 
scopic number of nodes coexists with finite components. 
The transition between these two phases occurs when the 



o o o 
o 




FIG. 1: Illustration of the linking process. Empty circles 
indicate active nodes, and filled circles depict inactive nodes. 
The initial state is a set of isolated nodes (left), while the 
final state is a regular random graph (right). As links are 
added, connected components form, and eventually, a single 
connected component spans the entire system. 



number of links exceeds a threshold. For example, when 
d = 3, the critical link density, Lg, is 



0.577200. 



(1) 



When d = 1 the final state consists of dimers, while when 
d = 2 the final state includes multiple rings [sol . [3l| . 
In both cases, the system does not condense into a sin- 
gle component. We study the more interesting situation 
d > 3 where the system undergoes a percolation transi- 
tion and the final state is a fully connected graph. 

As a result of the linking process, connected compo- 
nents undergo an aggregation process. In the traditional 
binary aggregation framework (3^ [ssj , connected com- 
ponents are characterized by the total number of nodes 
[22l424l |. In our process, however, this minimal descrip- 
tion does not lead to closed evolution equations. Instead, 
we describe clusters by the total number of nodes of de- 
gree 0, 1, . . ., d. This complete description naturally 
lends itself to closed evolution equations. We use the 
elegant Hamilton-Jacobi method to obtain a formal so- 
lution for the density of connected components, and find 
the percolation thresholds from this solution. 

We also use the degree distribution to illuminate the 
behavior during the late stages of the evolution. The 
regular random graph forms via a two-step process. 
First, the giant component takes over the entire sys- 
tem. This step occurs after T linking attempts, with 
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T ^ N / {hi NY~^. At this point, only a small minority 
of the nodes in the giant component remain active. The 
total number of active nodes, M , scales logarithmically 
with system size, 



M - (IniV) 



(2) 



Eventually all nodes become inactive, and the regular 
random graph is complete. The number of linking at- 
tempts required to reach this final state is proportional 
to the system size. Hence, the evolution of the regular 
random graph is much slower than the evolution of the 
classical random graph, where T is logarithmic in N [13] . 

The rest of this paper is organized as follows. The de- 
gree distribution is derived in section |lll We then show 
that the linking process is equivalent to a multivariate 
aggregation process (Sec. IIIip . A formal solution is ob- 
tained in Sec. II VI using the Hamilton- Jacobi method. The 
emergence of the giant component is discussed in section 
IVl In particular, the mass of the giant component is 
obtained as a function of time. The late-time kinetics 
including the emergence of a fully connected graph and a 
perfect regular random graph are described in section lVTl 
We discuss the results in section IVlIl In Appendix A, we 
rederive relevant properties of classical random graphs 
using the Hamilton- Jacobi method. 



The density of active nodes controls the process: the de- 
gree distribution obeys the rate equations 



drij 
~dt 




j < d, 
3=d. 



The initial condition is 



(5) 



(6) 



and the "boundary" condition n_i = ensures that 
dno/dt — —vriQ. Equations ([5|) are nonlinear because (j3]) 
is a two-body process. Yet, since the density u merely 
sets the overall linking rate, we can linearize Eqs. ([5]) by 
introducing the time variable 



dt' v{t') 



(7) 



or equivalently, dr /dt = v. In terms of the time variable 
T, the rate equations become 



drij 



nd- 



j < d, 
J^d. 



(8) 



We use the initial condition ^ and solve ([5]) recur- 
sively to find no = e^"^, ni — re^'^, n2 — ^r^e"'^, etc. 
Therefore, the degree distribution is a truncated Poisson 
distribution [SJI 



II. DEGREE DISTRIBUTION 

We study a random process that generates a regular 
random graph. The process starts at time t = with 
N disconnected nodes. At each elementary step, two 
nodes are chosen at random. If the degrees of the two 
nodes are both smaller than d, a new link connects them; 
otherwise, no link is added. Time is augmented by 2/N 
after each linking attempt, t ^ t + 2/N, so that every 
node participates in one linking attempt per unit time, on 
average. Linking is repeatedly attempted until a regular 
random graph forms (Fig. [T|). 

There are two types of nodes — active nodes with de- 
gree smaller than d and inactive nodes with degree equal 
to d. A new link between two active nodes of degrees 
i < d and j < d augments the degrees: 



(z,j) ^ (^ + l,J + l) 



(3) 



This random process occurs with unit rate. The linking 
process ([3]) transforms the system from a set of active 
nodes into a set of inactive nodes, as illustrated on Fig.[T] 
Let nj{t) be the degree distribution, that is, the frac- 
tion of nodes with degree j at time t. The degree dis- 
tribution is properly normalized, X)j=o "-i ~ ^^'^ 
partial sum gives the total density of active nodes, 



(9) 



for j < d. The density of inactive nodes, nd — I ~ ly, fol- 
lows from the normalization condition. From the defini- 
tions (dJ and ([7]) , the time variable r obeys the nonlinear 
ordinary differential equation 



dr 



(10) 



and 



with the initial condition r(0) = 0. Equations 
(fTO|l specify the degree distribution. 

The time evolution of the degree distribution, obtained 
by numerical integration of ([10]) for the case d = 3, is 
shown in Fig. [21 As expected, the density of isolated 
nodes, ?T.o(i), declines monotonically, while the density of 
inactive nodes, 71.3 (i), increases monotonically. Isolated 
nodes initially dominate, then nodes with degree one are 
most numerous, and eventually, inactive nodes take the 
lead for good. 

The average degree, (j) = X^j^o-^^i' characterizes the 
overall connectivity. Since every link connects two nodes, 
the density of links, L, is proportional to the average 
degree, L = (j)/2. By multiplying ([8]) by j and summing 
over j, we see that the average degree satisfies 



d^ 
dr 



(11) 



(4) Using the initial and final values, (j)o = and (j) 
we have the integral identity /„ dr v{t) = d. 
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FIG. 2: The degree distribution nj{t) versus time t for d = 3. 



III. MULTIVARIATE AGGREGATION 

Every node belongs to a connected component, or 
"cluster". Initially, there are N clusters, and eventu- 
ally, a single cluster remains. The number of clusters 
dwindles as a result of a cluster-cluster aggregation 35 1 
process: any link between two active nodes in separate 
clusters merges the two clusters. Since the aggregation 
rate depends on the total number of active nodes, char- 
acterization of clusters by their total size is insufficient. 
Adding the number of active nodes as a second variable 
does not lead to a closed description either, because the 
number of active nodes may change once a link is added. 
Therefore, we use a complete description and character- 
ize clusters by (d -I- 1)— dimensional vectors 



k = (fco, fci, • • • ,kd), 



(12) 



where kj is the number of nodes of degree j. The clus- 
ter size, k, is simply k = \k\ = k^ + ki + • • ■ + kd- Finite 
clusters are typically trees [s^ . Since there are k — 1 
links in a tree of size fc, we have the constraint 



2{k -l) = ki+ 2k2 



dkd 



In particular, minimal clusters with fc=(l,0, ••• , 0) 
trivially satisfy this condition. 

As a result of the linking process ([3]) , clusters undergo 
the multivariate aggregation process 



7 I 

I , m — > I + m + Bi + Bj, 



(13) 



with i < d and j < d. Here, the vector Bj describes the 
degree augmentation J — >■ j -t- 1, 

eo - (-1,1,0,0,...,0), 
ei = (0,-1, 1,0,. ..,0), 



The aggregation rate li x rrij in ([T3l) equals the product 
of the number of nodes of degree i in the first cluster and 
the number of nodes of degree j in the second cluster. 

The density c{k,t) of clusters with "size" k at time t 
obeys the rate equation 



dc{k) 
dt 



= '^'^hmjc{l)c{m)Sk,. 



i,j<d 




where v is the density of active nodes (|4]). Henceforth, 
the time dependence is implicit. The gain term directly 
reflects the aggregation process (|13|) . and the factor 1/2 
prevents over-counting. The loss rate is proportional to 
the total number of active nodes in the cluster, as well 
as the density of active nodes in the entire system. 

A crucial feature of the above rate equation is that the 
gain term and the loss term are different in nature. The 
gain term, which describes creation of clusters with finite 
size, is proportional to the product of cluster densities. 
The loss term, which describes loss of finite clusters due 
to linking, is proportional to the product of the cluster 
density and the total density of active nodes in the entire 
system, which does take into account the giant compo- 
nent, if one exists. For this reason, the rate equation 
is valid in the non-percolating phase as well as in the 
percolating phase. 

In terms of the time variable r, the loss term becomes 
linear 



dc{k) 



1 



limjc{l)c{m)6k. 



\kA c{k). 



dr 2iy^' 

i,j<d 

(14) 

The initial condition is c(fc,0) — i5fc,(i,o,...,o)- The rate 
equations for the cluster size distribution are nonlin- 
ear, in contrast with the linear equations ^ for the de- 
gree distribution. Moreover, the rate equations (|T4l) are 
substantially more challenging than those previously an- 
alyzed in studies of multivariate aggregation processes 
[37, 38]. 

The generating function approach is generally useful 
in aggregation problems [l^, S, IH, , and in our case, 
we must use the multivariate generating function 



C(a;,r)=^c(fc,T)a;^ (15) 

k 

Here x = {xq, xi, ■ ■ ■ , Xd) and x'' = Xq°Xi^ ■ ■ ■ x^"^ is a 
shorthand for the product of monomials. Multiplying 
(1141) by x'' and performing the summation over k, we 
find that the generating function evolves according to 




ed-i 



(0,0,. ..,0,-1,1). 



\j=o 

The initial condition is 

C{x, 0) — xo . 



d-l 
3=0 



dC 
' dxj 



(16) 



(17) 
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IV. HAMILTON-JACOBI FORMULATION 



for j < d, with function u defined by 



The generating function satisfies a non-linear Partial 
Differential Equation (PDE) with non-constant coeffi- 
cients. Crucially, equation is a first-order PDE, and 
hence, it can be reduced to a set of coupled Ordinary 
Differential Equations (ODEs) by applying the method 
of characteristics (ill, [4^] ■ This is a significant simplifi- 
cation as a PDE is essentially a collection of infinitely 
many ODEs. Moreover, in the present case, we conve- 
niently treat the governing equation as a Hamilton- Jacobi 
equation, thereby bypassing the method of characteris- 
tics [41I . |43| . Then, the canonical Hamilton equations 
coincide with the equations for the characteristics. 

First, we convert (jl6p into a Hamilton- Jacobi equation 



dC{x,T) 



H{x,VC,t) = 



(18) 



and view C{x,t) and H{x,VC,t) as "action" and 
"Hamiltonian" , respectively. Second, we treat x as "co- 
ordinate", and VC as "momentum", p — VC; in compo- 



nents, pj — for all j. The Hamiltonian is 



H{x,p,t) = ^XjPj 

j=o 



Mr) ' 



(19) 



where for convenience, we introduced a set of auxiliary 
functions. 



Pi 



-3 ' 



(20) 



for all j. We stress that the Hamiltonian (|19p is not 
a conserved quantity because of the explicit dependence 
on time t, which enters via v{t). The Hamilton- Jacobi 
equation (1181) is equivalent to the canonical Hamilton 
equations 



dxj 
d7 



dH 

dpj 



dr 



dH 



(21) 



Hence, the non-linear partial differential equation (|16p 
reduces to 2(d-|-l) coupled ordinary differential equations 
for the coordinates Xj and the momenta pj . 

From ([21]) and ([T9|) . the evolution equations for the 
momenta are 




Pj 



j < d, 
j^d. 



(22) 



The initial condition pj{0) — Sjfl follows from p = VC 
and (jl7p . and the boundary condition is p_i = 0. The 
structure of the linear equations ([22)) for pj is identical to 
that of equations ([5]) for nj , and the initial conditions are 



the same, too. By integrating equations 
we find 



P3 



recursively. 



(23) 



ni(r') 



nr') 



(24) 



or alternatively, du/dr = IVi/v. To determine pd^ we 
rewrite the last equation in (|22p as dpd/du = Pd~i, and 
therefore, pd = /g du' pd-i{u'). Formally, the momenta 
Pj resemble the degree distribution Uj. 

Using ([2T|) and (fT9)) . the evolution equations for the 
coordinates are 



dxj 
d7 



Vj+l 







j < d, 
J^d. 



(25) 



The initial condition is x{t = 0) = y. One of the co- 
ordinates does not change, Xd = yd- We may integrate 
(PS)) and solve for the final coordinate a; as a function of 
the initial coordinate y. Yet, as shown below, we need 
to solve the inverse problem — find y as a function of x. 

Fortunately, the auxiliary functions defined in (pUf en- 
able us to accomplish this task. We first note that at time 
T — 0, these sums coincide with the initial coordinates, 
Hj(r = 0) = Tjj. Using the evolution equations (1^^ and 
(I^SI) , we calculate the evolution equations for Hj , 



i-=3 

E {x, - ^x,+^ p,-j + iirP^-j-i - P^-j) 

-Xd Pd-j, 



1=3 



when j > 0. Here, the overdot is shorthand for time 
derivative, ' = A similar calculation shows that Hq is 
conserved. Ho = 0. Hence, the functions Hj obey 



dr 







-Xd Pd-j 



J > 0. 



(26) 



The initial condition is Hj(0) — yj. Since the coordinate 
Xd is a constant, integration of these evolution equations 
is immediate, and using the initial condition Hj(0) = yj, 
we find H 



yj = 



Ed 



J2t=j XtPt-j + Xd /J" dr'pd-jir') j > 0. 



(27) 



(The lower expression holds up to j = d] indeed, using 
Po = e~'^ we recover yd — Xd-) Equations (|27p . together 
with the known momenta (|23p , specify the initial coordi- 
nates in terms of the final coordinates. Remarkably, we 
formally solved the inverse problem! 

To complete the solution of the Hamilton- Jacobi equa- 
tions, we must specify the function u in (|23p . Using the 
definition ([M|) . we calculate the second derivative of u, 



d'^u _ d Ui 
dr^ dr v 



-Xd- 



Pd-i nd-i du 
V dr 



5 



Here, we used dv/dr = —Ud-i, as follows from ([8]). By 
multiplying the above equation by v, and by substituting 
dU and ([23| , we find that u satisfies a dosed second-order 
ordinary differential equation, 




du 



,d-l 



Xd- 



{d-l)\dT 



0. (28) 



The initial conditions are uifi) — and du(0)/dT — yi. 
This equation is linear if and only if c? < 2. 

However, we still do not know yi because the initial co- 
ordinates are coupled to the function u, according to (P7|. 
To find yi, we have to integrate equation p8|) . starting 
with trial initial conditions u(0) = and du{0)/dT = yl, 
and find the root of the integral equation 



Vi 



(29) 



The second order ODE ([25]) and the integral equation 
P9p specify the variable u and hence, the momenta pj 
and the initial coordinates yj via the explicit solutions 
and (|27p . We stress that the function u = u{x,t) 
depends on the coordinate x. 

Finally, we obtain the generating function (jlSp by eval- 
uating its complete derivative with respect to t. 



dC{x,T) 
d7 



dC dC dx 
dr dx dr 
lif V f du 
'2v ~ 2 \ dr 



-H + p 



dH 
dp 



(30) 



By integrating ([5(7| subject to the initial condition 
C{x, 0) — yo, and recalling the conservation law j/o — Hq, 
we get 



C{x,t) 



0=0 



dr' v{t') 



(31) 



This expression, supplemented by equations (|28p and 
((29| for u, as well as the explicit expressions (|23p for 
the momenta, constitutes a formal solution for the gen- 
erating function. Therefore, representation of the linking 
process as a multivariate aggregation process enables an- 
alytical treatment of the cluster-size distribution. 



V. GIANT COMPONENT 

By definition, the size distribution c(fe) describes finite 
connected components. However, the system may be in 
one of two phases — a non-percolating phase in which 
finite components contain all the mass, and a percolating 
phase in which a giant component, containing a fraction 
g of the entire mass, coexists with finite components. We 
expect that the giant component emerges suddenly, at a 



finite time tg (see Appendix A). 



The mass of the giant component g follows directly 
from the momenta when x — (1,1,..., 1). Indeed, using 
PS)) . we have 



E 

k 



kc{k) 



d 
3=0 



(32) 



or alternatively, g = I — yo, as follows from the conser- 
vation law Ho = yo- Hence, we substitute = 1 into 
equation (j28p . and obtain the evolution equation for u, 




^d-i 



du 



,d-l 



{d-l)\ dr {d-iy. 



^ 0. 



(33) 



The initial conditions are w(0) — 1 and du{0)/dT = j/i. 
Also, when x = (1,1,..., 1), the integral equation ([2^ 
becomes 



d-l 



(34) 



with Pj given by (j23p . 

First, there is the trivial solution u — t and yj — 1 for 
all j, for which the momenta equals the degree distribu- 
tion, Pj = rij. For this solution, g = 0, and the system is 
in the finite component phase. When d < 2, the trivial 
solution holds at all times. In the case d = 1 the final 
state is merely a collection of N/2 dimers; in the case 
d = 2, the final state includes only rings [33,[3l|. In both 
of these cases the system does not condense into a single 
component. 

When d > 3, there is a second nontrivial solution when 
T > Tg or equivalently, t > tg. To find the nontrivial solu- 
tion, we solved equation (I33p numerically, using Adams' 
method [i^ H^, starting with a trial initial condition 
m(0) = and du{0)/dT = yl, and then found the root 
of the integral equation ([M]) using the bisection method 
[47| . (Conveniently, the integral equations may be con- 
verted into first-order differential equations [4J].) We 
then obtained the mass of the giant component using 
([5^ . To evaluate time, t, and the average degree, (j), 
in terms of the modified time variable r, we simply inte- 
grated dt/dr — l/v and equation (fTTjl . 

For the case d — 3, we find the critical time 
Tg — 1.197395 that corresponds to 



tg = 1.243785. 



(35) 



When the giant component emerges, the average degree 
of a node is {i)g — 1.154399. Moreover, the critical den- 
sity of links quoted in ([l} is obtained from Lg = {j)g/2. 
Thus, the fraction fg of successful linking attempts is 



tn 



= 0.928135. 



(36) 



Generally, the percolation parameters characterizing the 
phase transition from the finite component phase to the 
giant component phase are nontrivial. 




FIG. 3: The giant component mass g{t) versus time t. Results 
of the Hamilton- Jacobi theory are shown for d = 3, 4, 5. Also 
shown as reference, is the solution to Eq. (|A7|) for the classical 
random graph, d = cxa 
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TABLE I: Percolation parameters for 3 < d < 7 and d — oo. 
Listed are the time, tg, the modified time, Tg, the average de- 
gree, {j)g = 2Lg, the fraction of successful linking attempts, 
fg , and the density of active nodes, Ug , when the giant com- 
ponent emerges. 

Table I lists the percolation parameters for 3 < d < 7. 
Remarkably, the restriction on the degree of a node has a 
small effect on the percolation parameters already when 
d = A. Moreover, as the maximal degree d increases, the 
percolation parameters quickly converge to the classical 
random graph values. 

Figure [3] shows the mass of the giant component for 
3 < < 5 and d — oo. At moderate times, the quantity 
g converges to the classical random graph behavior (jA7p . 
Yet, as will be shown in the next section, the long-time 
evolution of the regular random graph is much slower 
compared with the classical random graph. 

The function u at x — (1,1,...,!) also gives the total 
density of clusters, c, since c = C(l, 1, . . . , 1). The formal 
solution (|5T|l together with ([5^ give 

c(r) ^l-g{r)~\f^ dr' v{t') . (37) 

Figure |4] shows that the concentration is a smooth func- 
tion of time, even at the percolation point [4^. We 
also mention a useful relation between the cluster den- 
sity and the link density throughout the non-percolating 
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FIG. 4: The concentration c(t) versus time t for the case 
d = 3. 

phase. Every link reduces the number of clusters by one. 
Hence, dc/dt = —dL/dt, and consequently, c + L = 1 
when t <tg. In particular, Cg + Lg ^ I. 

Let Cfc be the density of clusters with total size fc, 
and let C(x) = J2k '^kX^ ^le the corresponding generating 
function. This single-variable generating function coin- 
cides with the multivariate counterpart (jl5l) when all the 
coordinates are identical 

e(x) = C(x,x, (38) 

Therefore, we have to vary x and obtain u = u{x) by 
solving the differential equation along with the in- 
tegral equation p9)) with Xi = x for all i. Finally, we 
substitute u into (|5T|) . We present results for the criti- 
cal behavior, Gg{x) = G{x,Tg), and focus on the large-A: 
behavior, that follows from (figure [S]) 

eg{x) ^cg-{l-x) + B{1 - x)3/2 + . . . , (39) 

as X ^ 1. The constant term is simply the critical con- 
centration, Cg — I — Lg = 1 — {.j)g/2 (Table I). The lin- 
ear term reflects the fact that finite clusters contain all 
the mass, C'(l) = J2k = 1. Finally, the leading sin- 
gular term (1 — x)^/^ shows that Ck ~ fc~^/^ as follows 
from the Taylor series, 

4v^^ r(fc-i-i) 

and the limiting behavior r{k + a)/T{k + b) ~ when 
fc — > (X). Hence, the critical cluster-size distribution has 
the power-law tail 

Ck^Ak-^^^, (40) 

with A — SB/Ay/n (see also Appendix A). Table II lists 
the prefactors i? for 3 < d < 7 and d = oo. As a critical 
phenomenon [4^, the evolving regular random graph is 
in the universality class of mean- field percolation [7| . 
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FIG. 5: The leading singular component of the generat- 
ing function Cg{x) at the critical point. The quantity 
A(a;) = Sgix) -Cg + {l~x) and the fit B{1 - a;)^/^ are plot- 
ted versus 1 — x. 
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0.962 


0.944 


0.943 


0.942809 



TABLE II: The prefactor B in ^ for 3 < d < 7 and d = oo. 
The last entry, B = 2\/2/3, is from Eq. (|A8| . 



VI. LARGE TIME BEHAVIOR 

We now turn to the large-time behavior. This behav- 
ior largely follows from the degree distribution We 
start with the crude estimate, t ~ Int, obtained by re- 
placing the right-hand-side of ((TU)) with the dominant 
exponential term, dr/dt = e~'^ . Then, we establish the 
corrections to this logarithmic behavior by keeping the 
leading term in (|10p 



'di 



{d-l)\ ■ 
From this equation, we find the relation 



r 



when t — >■ oo. Substituting this expression into © gives 
the leading asymptotic behavior of the degree distribu- 
tion 



(41) 



for j < d. The densities of active nodes with degree 
j < d ~ 1 decay algebraically with time, albeit with a 
logarithmic correction. Furthermore, the degree distri- 
bution is an increasing function of j. Since v ~ rid-i, 
the quantity v exhibits the universal power-law decay 



(42) 



The smaller the degree j, the faster the decay of rij. 
In particular, the density of isolated nodes exhibits the 
fastest decay. 



no 



(d- l)!t-i(lnt)-('^-i'. 



(43) 



Each isolated node represents a minimal cluster with size 
fc = 1. Therefore (|43|) gives the density of minimal 
clusters. The asymptotic behavior P5|) is much slower 
than the fast exponential decay found in ordinary ran- 
dom graphs [131 ■ Hence, the restriction on the degree 
leads to much slower long-time kinetics. 

The asymptotic results (02]) and (|^ provide use- 
ful information about the final stages of the evolu- 
tion. Amongst all clusters, minimal clusters survive the 
longest. Hence, when all minimal clusters disappear, a 
single cluster remains, and the giant component takes 
over the entire system. The time to reach a single com- 
ponent, T, grows with N according to 



N 



(In TV) 



d-i ' 



(44) 



as follows from the "depletion" criterion Nuq ^ 1. This 
time is not self-averaging: it fluctuates from realization 
to realization. In particular, the average (T) does not 
give the second moment (T^); namely, 



lim 



(r2 



> 1. 



(45) 



This breakdown of self averaging is in sharp contrast with 
the behavior of ordinary random graphs where T is a 
self-averaging quantity. The conclusion pS)) is based on 
the observation that a slight variation in the depletion 
criterion, say Nuq ~ 2, gives a different estimate for T. 

When the giant component contains all N nodes, the 
vast majority of the nodes are inactive. Yet, a small ma- 
jority of nodes remain active. Let Mj be the average 
number of nodes with degree j > when the graph be- 
comes fully connected for the first time [sOl- This quan- 
tity grows logarithmically with system size (figure [6]) 



Mj - (In NY, 



(46) 



with < j < d, as follows from Mj ~ Nnj{T), along 
with equations (j4T|) and (|44| . The majority of active 
nodes have the largest possible degree, j = d — 1. As 
stated in Eq. ([2]), the total number of active nodes, when 
the giant component first takes over the entire system, 
M, is logarithmic. 

To test these predictions, we performed Monte-Carlo 
simulations. In the simulations, we start with N isolated 
nodes. At each step, we pick two active nodes and con- 
nect them. Keeping track of the active nodes maximizes 
computational efficiency, and allows us to simulate large 
systems with N — 10^. We performed simulations for the 
case d — 3, and measured the average number of nodes 
with degree I and degree 2 when the graph condenses into 
a single component for the first time (figure [B]). Also, we 
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lo' 10^ 10^ lO'* 10^ 10^ 10^ 

N 



the leading asymptotic behavior (|48|) can be estabhshed 
recursively from equation (j49p . 

In general, we must classify clusters by their topol- 
ogy, and construct specific rate equations for each topol- 
ogy. The resulting rate equations have the same overall 
structure as (j49| and as a consequence, the asymptotic 
behavior appears to hold also when k > d. The pro- 
portionality constants depend on cluster topology, how- 
ever. Finally, we note that the criterion Nck ~ 1 shows 
that the asymptotic behavior is relevant only in a 
limited range of sizes, k with the logarithmic cut- 

off k^: :^ id — IniV. The cluster size distribution is 
sharply suppressed beyond the cutoff /c* . 

VII. DISCUSSION 



FIG. 6: The average number of active nodes when the graph 
becomes fully connected for the first time. Shown are Mi and 

1/2 

as a function of system size N , for the case d = 3. The 
results are from an average over lO'' independent simulation 
runs. 

verified that the time T fluctuates from realization to 
realization by measuring the ratio in (|45|) . 

The asymptotic behavior (PT|) indicates that there are 
several intermediate steps until the regular random graph 
fully forms. Isolated nodes disappear first, and at this 
time, there are (InA^)'^"^ active nodes. Nodes with de- 
gree 1 disappear next, and at this point, (InA^)''"^ ac- 
tive nodes remain. There are d such steps, and in the 
last step, nodes with degree d — 1 disappear. Now, the 
regular random graph is complete. The time it takes to 
form the regular random graph is linear in system size, 

tf ~ N, (47) 

as follows from the depletion criterion Nu ^ 1, and the 
algebraic decay Like the quantity T, the forma- 

tion time tf is a fluctuating quantity, even in the limit 
iV ^ oo. 

For completeness, we also mention the asymptotic be- 
havior of the cluster-size distribution, 

Cfc~t-i(lnt)-^-('^-i), (48) 

in the long-time limit (cZ = 2 is a special case of a model 
considered in fEU). For minimal clusters, ci = tiq, and 
this behavior follows from the density of isolated nodes. 
The behavior follows from the master equation 

^ = i ^ ImciCm-vkck, (49) 

l+tn—k 

for k < d. Except for the loss rate, which is proportional 
to the overall density of active nodes, i^, this rate equation 
is the same as (jAip . Clusters with k < d contain only 
active nodes, and hence, the loss term is proportional 
to cluster size, k. Starting with ci ~ i~^(lnt)~(''~^\ 



In conclusion, we studied the evolution of random 
graphs with a bounded degree using a dynamic link- 
ing process. First, a giant component, which contains 
a macroscopic number of nodes, emerges in a finite time. 
We obtained the nontrivial thresholds for this percola- 
tion transition. Then, the giant component takes over 
the entire system. The evolution still proceeds as some 
nodes have yet to reach the maximal degree. While the 
number of linking attempts before the giant component 
emerges is a self-averaging quantity, the number of steps 
it takes to form a fully connected graph or to complete 
the regular random graph are both fluctuating (non-self- 
averaging) quantities. 

The regular random graph is completed via multiple 
steps. In the first step, nodes with degree disappear 
and at the point, the graph is fully connected. Then, 
nodes with degree 1 disappear, etc. The total number of 
nodes with degree d — 1 diminishes by a factor In N in 
each step, 

{InNf-^ -> (InA^)''-^ ^ ^ (InA^) -> 0. 

Similarly, the number of nodes of degree d—2 shrinks ac- 
cording to (IniV)*^"^ — ^ (IniV)''"^ ■ ■ ■ . Furthermore, 
the asymptotic behavior (l49l) indicates a multitude of 
time scales. Clusters of size k disappear at time 

N 

^ (IniV)Md-i)' 

for all k ^ k^,. The time T = Ti is simply the largest in 
a series of time scales, Ti > T2 > > ■ ■ ■ . Therefore, 
multiple finite-size scaling laws characterize kinetics of 
regular random graphs. 

The evolution equation for the cluster-size density is 
not a closed equation. To address this challenge, we 
introduced the multivariate size distribution, where the 
number of nodes with a given degree is specified. For 
this quantity, the evolution equations are closed. Us- 
ing the Hamilton- Jacobi method, also useful for analyz- 
ing large fiuctuations in population dynamics [52| - [55j , we 
constructed a formal solution for the multivariate size 
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distribution. Moreover, the cumbersome rate equations 
reduce to a single second order differential equation, and 
numeric integration of this equations gives the percola- 
tion thresholds with, essentially, arbitrary accuracy. 

Our analysis shows that the rate equation description 
extends to a broader class of evolving random graphs. 
Certainly, the multivariate aggregation analysis can be 
generalized to situations where the linking rate is degree 



dependent. The multivariate aggregation framework can 
also be used to study the structure of clusters and 
in particular, cycles. 

We thank Wolfgang Losert for useful discussions. 
This research was supported by DOE grant DE-AC52- 
06NA25396 and NSF grant CCF-0829541. 
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Appendix A: Classical Random Graphs 

This appendix briefly describes how to obtain relevant 
properties of the classical random graph, d = oo, us- 
ing the Hamilton- Jacobi formalism. The density of con- 
nected components with size k obeys [l^l 



dck _ 1 
'dt ^ ~ 



kch 



(Al) 



l-^m—k 



with the initial condition Cfe(O) = (5^,0 . The generating 
function Q{x,t) — '^i^ Ck{t)x^ satisfies 



96 o»e 1 / ae 



dt^"" dx 2 r dx 



(A2) 



with the initial condition C(a;,0) — x. 



Using the momentum P — the Hamiltonian is 



dx 
1 



H = xp- ^{xpf. 



(A3) 



The Hamiltonian does not depend explicitly on time. 
Hence, H is a conserved quantity and consequently, xp 
is also constant. The evolution equations are 

^ ^ x{l- xp), ^ ^ -p{l- xp), (A4) 

with the initial conditions a;(0) = y and p{0) — 1. We 
can verify that d{xp)/dt = 0, and therefore, xp = y. 
By integrating (IA4p . we obtain the coordinate and the 
momentum. 



x^ye^^-y^', p^e-^'^-y^'. 



The evolution equation for generating function is 



(A5) 



Hence, the generating function is a quadratic function of 
the initial coordinate. 



eix,t) 



1 2 

y-^y t, 



(A6) 



=(i-y)t 



with X = y e'- 

The mass of the giant component, g = 1 — X^fc^'^fcj 
follows from the momentum p at x = 1. The conservation 
law xp — y implies that p = y if a: = 1. In analogy with 
(|32l) . we have g = 1 — y, and using (jA4p . we find a closed 
equation. 



1 



-gt 



(A7) 



There is a trivial solution g — 0, valid at all time, and a 
second nontrivial solution, g > 0, valid only when t > 1. 
The latter is physical when t > 1, and hence, the giant 
component emerges at time tg = 1. 

To obtain the critical size distribution, that is, 
the size density at t — tg, we note that the crit- 
ical generating function is given by Gg{x) ~ y — ^y^ 



^ We then focus on the be- 

1. We substitute the expansion 



with X — ye 
havior for x — 

1 — y = ai(l — x)^/^ + 02(1 — y) + • • • into the equation 
X = ye^~y , and solve for the coefficients ai and 02. This 
calculation gives 



l-y=[2(l-x)]i/2__(i_^) + ... . 

This expansion, together with Eq. (|A6[) . gives the leading 
behavior of the critical generating function. 



Q,{x)^\-{l-x)+^-^{l-x)^'^ + 



(A8) 



in the limit x — 1. Thus, the critical size distribution 
has the power-law tail [STj 



Cfe 



/2^ 



(A9) 



as follows from the leading singular component (1 — o;)'^/^. 



